module globals
	use myfuncs
    use dotops
	use ieee_arithmetic
    
	character(len=80)	:: folder="c:/dropbox/mystuff/heavymean/2020/swind/"

! comment in the (k,n0) block that we want 
	
	integer, parameter	:: k=8
	real				:: swsmall=0.5, swflat=1.8
	integer, parameter	:: n0=50
	real, parameter		:: varestfac=2.0

	!integer, parameter	:: k=4
	!real				:: swsmall=0.5, swflat=1.5
	!integer, parameter	:: n0=50
	!real, parameter		:: varestfac=2.0

	!integer, parameter	:: k=4
	!real				:: swsmall=0.5, swflat=1.5
	!integer,parameter	:: n0=25
	!real, parameter		:: varestfac=2.0

	!integer, parameter	:: k=12
	!real				:: swsmall=0.7, swflat=2.5
	!integer, parameter	:: n0=50
	!real, parameter		:: varestfac=2.0

	real				:: maxconsttail=1E-10
	real, parameter		:: cvglobal=1.0
	logical, parameter	:: checkflag=.false.
	logical, parameter	:: tstconst=.true.			! if true, constrain tests to ensure coherence of CIs and p-values

	integer, parameter	:: nth=6, nx0=6, nth_one=3, nx0_one=3, nth_two=nth
	real, parameter		:: xsimin=-0.4 ,xsimax=.5
	real, parameter		:: xsimin1=-0.5, xsimax1=0.5
	real, parameter		:: impvar=2.0
	real				:: level=0.01, levelfudge=0.001
	
	integer, parameter	:: npb=128, nblocks=5000, nbi=(npb-1)/64+1, nfpc=npb/nbi
	integer, parameter	:: nsim=npb*nblocks
	logical, parameter	:: lowramflag=(k>12)

	integer				:: stage=0
	
	real				:: Y0s(0:k,npb,nblocks)
	real				:: densimp(npb,nblocks) , falt(npb,nblocks)

	integer(8)			:: credflags(nbi,npb,nblocks),ctest(nbi,npb,nblocks)
	integer				:: pairbind(nblocks)				! starting index minus one in pairs list for each block
	
	integer				:: nimp
	real, allocatable	:: thlistimp(:,:)
	
	integer				:: n_sample
	
	type ptype			! stats for pairs
		real	:: f1s,f1b,Zsum,Z(k,2)
	end type
	
	type(ptype), allocatable	:: pairsdat(:)
	integer, allocatable		:: pairsind(:,:)
	
	integer, parameter	:: n0max=2000
	integer				:: nsimeff, n0active, n0active_stage1
	real, allocatable	:: lam(:), lam_stage1(:), thlist0_stage1(:,:)
	real				:: thlist0(nth,n0max)
	
	integer, parameter	:: nGQ=40
	real				:: GQxw(nGQ,2)
	
end module

module basicmodule
	use globals
	implicit none

	contains
	
! ************************************************************************
! ***********  Generic Routines ******************************************	
! ************************************************************************
	
	
	function expecval(n,xsi) result(val)
        integer :: n
        real    :: xsi,val
        val=(gamma(n - xsi)/Gamma(real(n)) - 1)/xsi
	end function

	function getej(n,msx) result(val)
		integer	:: n
		real	:: msx(3),val
		val=msx(2)*(expecval(n,msx(3))+msx(1))
	end function

	function getdens_base(msx,Y) result(val)
        real    :: y(k),msx(3),m,sig,xsi,val
        real    :: zk
        m=msx(1); sig=msx(2); xsi=msx(3)
		if(m==99) then
			print *,"error: tiny tail getdens_base evaluation"
			val=0
			stop
			return
		endif
		if(Y(1)==1E-10) then
			val=0
			return
		endif
        zk=1+xsi*(y(k)/sig-m)
        if(zk<0) then
            val=0
            return
		endif
        if(1+xsi*(y(1)/sig-m)<0) then
            val=0
		else
			val=exp(-k*log(sig)-zk**(-1/xsi)-(1+1/xsi)*log(product(1+xsi*(y/sig-m))))	
		endif
		if(.not. isfinite(val)) then
			val=0
		endif
	end function
	
	
	function getmeanadj(msx,Yk) result(val)
        real    :: Yk,msx(3),m,sig,xsi,zk,val
		if(msx(1)==99 .or. Yk==maxconsttail) then
			val=0
			return
		endif
        m=msx(1); sig=msx(2); xsi=msx(3)
        zk=1+xsi*(yk/sig-m)
        if(zk<0) then
           val=0
        else
            val=(sig/xsi)*zk**(-1/xsi)*(zk/(1-xsi)+xsi*m-1)
        endif
		if(.not. isfinite(val)) val=0
	end function
	
	
	function getdens0(msx,Y) result(val)
        real    :: y(0:k),msx(3),val
		real	:: mv(2)
		if(msx(1)==99) then
			if(Y(1)==Y(2)) then
				val=exp(-0.5*(y(0)+k*maxconsttail)**2/0.5)/sqrt(0.5)
			else
				val=0
			endif
			return
		endif
		if(Y(1)==Y(2)) then
			val=0
			return
		endif
		mv=[getmeanadj(msx,Y(k)),0.5]
		val=getdens_base(msx,Y(1:))*exp(-0.5*(y(0)+mv(1))**2/mv(2))/sqrt(mv(2))
	end function

	function getdensimp(msx,Y) result(val)
        real    :: y(0:k),msx(3),val,mv(2)
		if(msx(1)==99) then
			if(Y(1)==Y(2)) then
				val=exp(-0.5*(y(0)+k*maxconsttail)**2/impvar)/sqrt(impvar)
			else
				val=0
			endif
			return
		endif
		if(Y(1)==Y(2)) then
			val=0
			return
		endif
		mv=[getmeanadj(msx,Y(k)),impvar]
		val=getdens_base(msx,Y(1:))*exp(-0.5*(y(0)+mv(1))**2/mv(2))/sqrt(mv(2))
	end function
	

	function getdens0_stage1(msx,Y,Zsum,var) result(val)
        real    :: msx(3),Y(k),Zsum,var,val,ma
		ma=getmeanadj(msx,Y(k))
		val=getdens_base(msx,Y)*exp(-0.5*(Zsum+ma)**2/var)/sqrt(var)
	end function
		
	function getdens0p_stage2(th,p) result(val)
		type(ptype)	:: p
        real    :: th(6),val,ma
		ma=getmeanadj(th(1:3),p%Z(k,1))-getmeanadj(th(4:6),p%Z(k,2))
		val=getdens_base(th(1:3),p%Z(:,1))*getdens_base(th(4:6),p%Z(:,2))*exp(-0.5*(p%Zsum+ma)**2)
		ma=getmeanadj(th(4:6),p%Z(k,1))-getmeanadj(th(1:3),p%Z(k,2))
		val=val+getdens_base(th(4:6),p%Z(:,1))*getdens_base(th(1:3),p%Z(:,2))*exp(-0.5*(p%Zsum+ma)**2)
		val=0.5*val
	end function

	function getdens0p(th,p) result(val)
		type(ptype)	:: p,pc
        real    :: th(6),val,var,Zsum,om
		integer	:: j
		if(stage==1) then
			Zsum=p%Zsum-sum(p%Z(:,2))
			var=1+sum(p%Z(:,2)**2)			
			val=getdens0_stage1(th(1:3),p%Z(:,1),Zsum,var)
			val=val/p%f1s
		else
			val=getdens0p_stage2(th,p)
			val=val/p%f1b
		endif
	end function

	subroutine setYfromeps(msx,eps,Y)
		real	:: msx(3),eps(0:k),Y(0:k),z(1:k)
		integer	:: i
		if(msx(3)==99) then
			Y(0)=sqrt(impvar)*eps(0)-k*maxconsttail
			Y(1:k)=maxconsttail
			return
		endif
		
  		z(1)=eps(1)
		do i=2,k
			z(i)=eps(i)+z(i-1)
		enddo
		z=z**(-msx(3))
		z=(z-1)/msx(3)
		Y(1:k)=msx(2)*(z+msx(1))
		Y(0)=sqrt(impvar)*eps(0)-getmeanadj(msx,Y(k))
	end subroutine
	
end module


	